home *** CD-ROM | disk | FTP | other *** search
/ Nautilus 1992 July / Nautilus-3-8 / Nautilus-3-8.bin / Tools & Utilities / Techy Stuff / Source ƒ / sky source ƒ / JUP.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  3.8 KB  |  194 lines
  1. #include "sky.h"
  2.  
  3. extern struct juptab
  4. {
  5.     float f[6];
  6.     char c[3];
  7. } juptab[];
  8.  
  9. jup()
  10. {
  11.     double pturbl, pturbb, pturbr;
  12.     double lograd;
  13.     double dele, enom, vnom, nd, sl;
  14.     double q0, v0, t0, m0, j0 , s0, u0, n0;
  15.     double lsun, elong, ci, dlong;
  16.     double planp[9];
  17.     struct juptab *pp = &juptab[0];
  18.     double olong;
  19.     double grin;
  20.     double temp;
  21.  
  22. /*
  23.  *    The arguments nnd coefficients are taken from
  24.  *
  25.  *    Here are the mean orbital elements.
  26.  */
  27.  
  28.     object = "Jupiter     ";
  29.     ecc = 0.04833748 + 1.63903e-4*capt - 0.4642e-6*capt2 - 1.593e-9*capt3;
  30.     incl = 1.3086429 - 0.005696*capt + 3.89e-6*capt2;
  31.     node = 99.4431901 + 1.0105300*capt + 3.522e-4*capt2 - 8.51e-6*capt3;
  32.     argp = 13.4119487 + 0.21344495*capt + 7.466e-4*capt2 - 3.7946e-6*capt3;
  33.     anom = 225.3350378 + 0.0830853474*eday - 8.332e-4*capt2 + 3.9876e-6*capt3;
  34.     motion = .083091212;
  35.     mrad = 5.202803;
  36.  
  37.     incl *= radian;
  38.     node *= radian;
  39.     argp *= radian;
  40.     anom = fmod(anom, 360.)*radian;
  41.     motion *= radian;
  42.  
  43. /*
  44.  *    Longitudes of perturbing planets,
  45.  *    They are of epoch Jan 0.5, 1900, and are
  46.  *    referred to the fixed qquinox of that date.
  47.  */
  48.  
  49.     q0 = 178.179 + 4.092338817*eday;
  50.     v0 = 342.767 + 1.602130491*eday;
  51.     t0 =  99.697 + 0.985609114*eday;
  52.     m0 = 293.748 + 0.524032950*eday;
  53.     j0 = 238.050 + 0.083091230*eday;
  54.     s0 = 266.280 + 0.033459743*eday;
  55.     u0 = 243.370 + 0.011731421*eday;
  56.     n0 =  85.183 + 0.005987356*eday;
  57.  
  58.     q0 *= radian;
  59.     v0 *= radian;
  60.     t0 *= radian;
  61.     m0 *= radian;
  62.     j0 *= radian;
  63.     s0 *= radian;
  64.     u0 *= radian;
  65.     n0 *= radian;
  66.  
  67.     grin = 5.*s0 - 2.*j0 - 8079.0*radsec*capt;
  68.  
  69.     planp[1] = q0;
  70.     planp[2] = v0;
  71.     planp[3] = t0;
  72.     planp[4] = m0;
  73.     planp[5] = j0;
  74.     planp[6] = s0;
  75.     planp[7] = u0;
  76.     planp[8] = n0;
  77.  
  78. /*
  79.  *    Computation of long period terms affecting the mean anomaly.
  80.  */
  81.  
  82.     anom += 0.
  83.     +(1192.85-6.076*capt-0.0400*capt2+0.00075*capt3)*radsec*sin(grin)
  84.     +(-23.80-0.192*capt+0.0226*capt2-0.00080*capt3)*radsec*cos(grin)
  85.     +(-11.04-0.060*capt-0.0072*capt2+0.00021*capt3)*radsec*sin(2.*grin)
  86.     +(1.44-0.086*capt+0.0004*capt2+0.00006*capt3)*radsec*cos(2.*grin)
  87.  
  88.      + (8.22-0.120*capt-0.002*capt2)*radsec*sin(2.*j0-6.*s0+3.*u0)
  89.      + (0.55 + 0.420*capt - 0.0079*capt2)*radsec*cos(2.*j0-6.*s0+3.*u0)
  90.     ;
  91.  
  92. /*
  93.  *    Computation of elliptic orbit.
  94.  */
  95.  
  96.     enom = anom + ecc*sin(anom);
  97.     do {
  98.         dele = (anom - enom + ecc * sin(enom)) /
  99.             (1. - ecc*cos(enom));
  100.         enom += dele;
  101.     } while(fabs(dele) > 1.e-8);
  102.     vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  103.             cos(enom/2.));
  104.     rad = mrad*(1. - ecc*cos(enom));
  105.  
  106. /*
  107.  *    Perturbations in longitude.
  108.  */
  109.  
  110.     pturbl = 0.
  111.         +(-83.79-1.222*capt+0.0097*capt2)*sin(grin-j0)
  112.         +(137.08-1.508*capt-0.0069*capt2)*cos(grin-j0)
  113.     ;
  114.  
  115.     for(;;){
  116.         if(pp->c[2]==0){
  117.             pp++;
  118.             break;
  119.         }
  120.         temp = planp[pp->c[2]]*pp->c[0] + j0*pp->c[1];
  121.         pturbl += (pp->f[0]+pp->f[1]*capt+pp->f[2]*capt2)*sin(temp)
  122.             + (pp->f[3]+pp->f[4]*capt+pp->f[5]*capt2)*cos(temp);
  123.         pp++;
  124.     }
  125.  
  126. /*
  127.  *    Perturbations in latitude.
  128.  */
  129.  
  130.     pturbb = 0.;
  131. /*
  132.     for(;;){
  133.         if(pp->f[0]==0.){
  134.             pp++;
  135.             break;
  136.         }
  137.         pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*j0 + pp->c[1]*planp[pp->c[2]]);
  138.         pp++;
  139.     }
  140. */
  141.  
  142. /*
  143.  *    Perturbations in log radius vector.
  144.  */
  145.  
  146.     pturbr = 0.;
  147. /*
  148.     for(;;){
  149.         if(pp->f[0]==0.){
  150.             pp++;
  151.             break;
  152.         }
  153.         pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*j0 + pp->c[1]*planp[pp->c[2]]);
  154.         pp++;
  155.     }
  156. */
  157.     pturbr *= 1.e-6;
  158.  
  159. /*
  160.  *    reduce to the ecliptic
  161.  */
  162.  
  163.     olong = vnom + argp + pturbl*radsec;
  164.     nd = olong - node;
  165.     lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
  166.  
  167.     sl = sin(incl)*sin(nd);
  168.     beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
  169.  
  170.     lograd = pturbr*2.30258509;
  171.     rad *= 1. + lograd;
  172.  
  173. /*
  174.  *    Compute motion for planetary aberration.
  175.  */
  176.  
  177.     temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
  178.     ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
  179.     bdot = temp*sin(incl)*cos(lambda-node);
  180.     rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
  181.  
  182. /*
  183.  *    Compute magnitude.
  184.  */
  185.  
  186.     mag = -8.93;
  187.  
  188.     semi = 98.57;
  189.  
  190.     helio();
  191.     geo();
  192.  
  193. }
  194.